if(1) 
  nLevel = 1;
  
  N = 64;
  ntotal = N*N;
  GRD = reshape(1:ntotal,N,N);

  display('Constructing Sample List first');
  tic

  SampleList = ConstructSampleList( N, nLevel );

  toc
  display(['Constructing Random Kernel upto Level ' int2str(nLevel)]);
  tic

  % First, construct a total random matrix, and then
  % substitute the low-rank part by low rank random matrix and
  % symmetrize it.
  Kernel = randn(ntotal, ntotal);

  nRealRank = 10;
  nWidth = N / 4;
  for iLevel = 1 : nLevel
    if( iLevel > 1 )
      nWidth = nWidth / 2;
    end
    for iSource = 1 : SampleList{iLevel}.nSource
      SourceBlock = SampleList{iLevel}.Source(:, iSource);
      Source = [(SourceBlock(1)-1)*nWidth+1, ...
	(SourceBlock(1))*nWidth; ...
	(SourceBlock(2)-1)*nWidth+1, ...
	(SourceBlock(2))*nWidth];
      for iTarget = 1 : SampleList{iLevel}.nTarget
	TargetBlock = SampleList{iLevel}.SourceTarget(:, iTarget, ...
	  iSource );
	Target = [(TargetBlock(1)-1)*nWidth+1, ...
	  (TargetBlock(1))*nWidth; ...
	  (TargetBlock(2)-1)*nWidth+1, ...
	  (TargetBlock(2))*nWidth];
	SourceInd = reshape( GRD(Source(1,1):Source(1,2), ...
	  Source(2,1):Source(2,2)), nWidth * nWidth, 1 );
	TargetInd = reshape( GRD(Target(1,1):Target(1,2), ...
	  Target(2,1):Target(2,2)), nWidth * nWidth, 1 );
	Kernel(SourceInd, TargetInd) = ...
	  randn(nWidth*nWidth,nRealRank) * ...
	  randn(nRealRank,nWidth*nWidth) + ...
	  1d-10 * randn(nWidth*nWidth, nWidth * nWidth );
      end
      end
    end

    Kernel = (Kernel + transpose(Kernel))/2;
  toc

end

if(1)
  display(['Constructing H matrix upto Level ' int2str(nLevel)]);
  tic

  MatDiag = ConstructHMatrix( Kernel, nLevel );
  toc

end


if(1)
  display(['Verifying the diagonal']);

  KernelDiag = reshape(diag(Kernel), N, N);
  err = norm( KernelDiag - MatDiag, 'fro' ) ./ norm( KernelDiag, 'fro' )
end
